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Introduction 

Scientific research often involves the measurement 
of two sets of variables to determine the functional re- 
lationship between the two variables. Because of the 
nature of experiments, the ordinate variable usually 
contains random errors. In those instances in which 
the underlying form of the relationship between the 
two variables is known, the relationship can be approx- 
imately obtained through the use of a mathematical 
tool such as the method of least squares (ref. 1). If the 
underlying relationship is unknown or if the relationship 
is known but is too complicated to be easily computed, 
the function is often approximated with polynomials 
and other simple functions. 

If the error in the ordinate data is negligible or 
zero, intermediate values can be calculated by inter- 
polating. Various spline interpolation functions have 
found widespread use over the years, including cubic 
splines (ref. 2), basis splines (B-splines) (ref. 3), and 
splines under tension (ref. 4). The last type of spline 
was devised to overcome the occasional occurrence of 
spurious local behavior of the cubic spline. The spline 
under tension reduces the deviant behavior by applying 
a constant tension to the cubic spline over the entire 
data range. The disadvantage to this approach is that 
the tension is also applied in those regions where the 
user would be satisfied with an untensioned spline. 

A lesser known tension spline developed by Spath 
(ref. 5), called the “rational spline,” allows the user to 
specify different tension values between each adjacent 
pair of knots. This is accomplished by introducing into 
the cubic terms of the spline a denominator containing a 
different parameter for each knot interval; depending on 
the parameter value, the rational spline can be made to 
behave like a cubic or a linear function. In this manner 
the behavior of the rational spline between each pair of 
knots can be adjusted by trial and error. Recently, Frost 
and Kinzel (ref. 6) introduced an automatic adjustment 
scheme which varies the tension parameter for each 
interval until the maximum deviation of the spline from 
the line joining the knots is less than or equal to a user- 
specified amount. This procedure frees the user from 
the drudgery of adjusting individual tension parameters 
while still giving him control over the local behavior of 
the spline. 

If the ordinates of the data contain errors, algo- 
rithms are available for finding least-squares approxi- 
mations by using cubic splines (refs, 7 and 8). Although 
the cubic-spline approximations provide smooth repre- 
sentations of the data, in some instances the approx- 
imations exhibit the same deviation that occurs with 
interpolating splines. For this reason, a spline approx- 
imation method which smooths the data without fluc- 
tuating wildly would be a very useful tool. 


This paper presents an algorithm for weighted least- 
squares approximations with rational splines. The al- 
gorithm contains the automatic adjustment procedure 
of Frost and Kinzel (ref. 6) for determining the tension 
parameters and a constrained weighted least-squares so- 
lution comparable to that found in reference 8 for es- 
timating spline coefficients. First, the derivation and 
description of the algorithm are presented. These are 
followed by an illustrative example comparing rational- 
spline and cubic-spline representations of the data. 

Derivation of Equations 

Let the set of n data points be represented by 
(x t ',y t ), where i = 1,2, ...,n and X\ < < . . . < 

x n , (A list of symbols used in this paper appears 
after the references.) Let Wi (for i = 1,2, . . . , n) be 
a set of positive weights indicating confidence in the 
corresponding ordinate values. Select / knots of the 
spline, the abscissas are Xfc, where k = 1,2,...,/ and 
Xi = x x < x 2 < . . . < xi-i <xi- x n . The ordinates of 
the knots are y k \ because the spline will pass through 
the knots, the knot ordinates are unknown at this point. 
Also, define dx k = x/c+i — x k for k = 1,2,...,/ — 1 
to be the length of each subinterval defined by two 
consecutive knot abscissas. The rational spline on 
interval k (k = 1, 2, ...,/ — 1) is defined in references 5 
and 6 to be 

ti 3 

F k (x) = A k u + B k t + C k p ^ t + 1 

+ Dk p kU+l ( x k< x < x k+ 1) (1) 

where 

P k tension parameter 

A ki B k ,C k , D k unknown coefficients 
u — (xfc+i - x)/dx k 

t = (x - x k )/dx k — l-u 

Equation (1) is defined for all independent variables 
x in the data range if the tension parameter P k is 
restricted to P k > -1. If P k is set to zero, equation (1) 
reduces to a cubic-spline function. As P k increases 
from zero, the cubic terms decrease and F k (x) tends 
to the equation of the line joining the knots at x k and 
Xfc+i. This characteristic of F k (x) is exploited in the 
automatic adjustment algorithm devised by Frost and 
Kinzel (ref. 6) to determine the value of P k . 

Evaluation of equation (1) for each subinterval re- 
quires knowledge of the four coefficients A ki B ky C k , 
and D k . This means that for / knots (equivalently, 
l — l subintervals), 4/ — 4 coefficients must be esti- 
mated by the method of least squares. The magni- 
tude of the estimation problem can be reduced by writ- 
ing F k (x) in terms of the unknown function values 



{yk for k = 1,2,...,/) and second derivatives (y'£ for 
k = 1, 2, ...,/ — 1) at the knots. To do this, first evalu- 
ate F k {x) at the knots x k and Xk+i and set respective 
results equal to y k and y k +x. If x — x k) then u = 1, 
t = 0, and 

Vk = F k (x k ) = A k + C k (2) 

If x = x k +i, then u = 0, t = 1, and 

Vk+i = F k (x k +i) = B k + D k (3) 

Next, differentiate F k {x) in equation (1) twice, evaluate 
the result at x k and x^+i, and set respective results 
equal to y f £ and y k+1 . The second derivative is 


^k u3 + 6P k {P k t + l)u 2 + 6 (P k t + l) 2 u 
kU k dx\{P k t + 1)3 

D 2P 2 / 3 + 6 P k (P k u + l)/ 2 + 6 (P k u + 1 )H 
dx 2 k (P k u+ l) 3 


Evaluating equation (4) at the knots yields 


(4) 


and 


yZ = FZ(x k ) = C k 


2 Pi +6 P k +6 

dx l 


(5) 


y'k+i = Fk(*k+ 1) = D k 


2 Pi + 6P fc + 6 

dx l 


(6) 


Solving equations (2), (3), (5), and (6) for A k , B ky C k , 
and D k yields the relationships 


where 


A k =Vk — Hky'k 

(7) 

Bk = 2/fc+i - HkVk+i 

(8) 

C k = H k jfi 

(9) 

Dk = H k y'k + 1 

(10) 

dx? 

k 2 (P fc 2 + 3 P k + 3) 

(11) 


Substituting equations (7) through (11) into equa- 
tion (1) and rearranging leads to the following rational- 
spline form: 


F k (x) = uy k + H k y’l + fy fc+1 

(12) 

Equation (12) shows the direct dependence of the ra- 
tional spline on the function and the second derivatives 
at the knots. If P k is set to zero, equation (12) re- 
duces to the equation of a cubic spline given in refer- 
ence 8. In this form, only the 21 parameters y k and y f £ 
for k = 1,2,...,/ need to be estimated from the data. 


Since equation (12) does not explicitly depend on 
the first derivatives at the knots, there is no assurance 
that the first derivative is continuous at the interior 
knots. To overcome this deficiency, the first derivatives 
must be constrained to be continuous at these knots. 
The first derivative of equation (12) is 



HjL 

dx k 


, gfc+l j Pk 
dx k dx k 


3u 2 ( P k t + 1) + u 3 P k 

( p k t + 1) 2 

~ 3t 2 (P fc u+l)+t 3 P fc 

(P k u + l ) 2 


The continuity constraints to be imposed are 


y'Ui 

(13) 


F 'k (*fc+i) = K + 1 (zfc+i) (* = 1, 2, .... Z - 2) (14) 


Using equation (13) for F^(x k + j) and F/. +1 (x fc+1 ) in 
equation (14) and rearranging yields the following Z - 2 
constraint equation (for k = 1, 2, ...,/ — 2): 


J£ 5i ' + (*s + 


J_Y 


This equation contains the six unknowns y and y " at 
the knots x k , x^+i, and x k + 2 - 

Equation (12) evaluated at the data points along 
with the constraint equation (15) form the basis for 
finding a rational-spline approximation. The approx- 
imation is found by solving a constrained weighted 
least-squares problem. In order to solve this problem, 
the following matrices are defined. Let Y be the 2/ 
column vector Y = (j/i,j/", fayZ, ■ . • , yi, y") T , where 
superscript T indicates matrix transpose. Let E be the 
n X 21 matrix containing cofactors of y k and y% (for 
k = 1,2,...,/) in equation (12) evaluated at the data 
abscissas x* (for i = 1, 2, . . . , n); each row corresponds 
to one data point. Because the coefficients [y k ,y k ) 
change from subinterval to sub interval, the matrix E 
has the following overlapping block structure: 



Details of the nonzero entries in E (those in the blocks 
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above) are given in appendix A. Let Y be the n X 1 
column vector containing the n data ordinates yi (for 
i = 1,2 , . . . , n); then the scalar equation yi = Fk{x{) 
for Xfc < < Xfc-f i and i = 1, 2, . . . , n can be written 

as the matrix equation 

Y = EY (16) 

Finally, let W be an n X n diagonal weight matrix 
containing the weights Wi (for i = 1 , 2, . . . , n) on the 
diagonal. 

The constraint equation (15) can be written in ma- 
trix form as 

SY — 0 (17) 

where S is the l — 2xl matrix of cofactors of y k and y f £ 
(k = 1, 2, ...,/) in equation (15). Because equation (15) 
relates y k and y % at three consecutive knots, S has the 
structure 


X 

X 

X 

X 

X 

X 

0 

0 

0 




O’ 

0 

0 

X 

X 

X 

X 

X 

X 

0 

0 

• 


0 

0 

0 

0 

0 

X 

X 

X 

X 

X 

X 

0 


• 

0 

0 

. 

. 

0 

X 

X 

X 

X 

X 

X 

0 

0 

0 

0 

• 


0 

0 

0 

X 

X 

X 

X 

X 

x _ 


The nonzero entries in S are defined in appendix A. 

With these matrix definitions the constrained least- 
squares problem requires the minimizing of (Y — 
EY) t W(Y —EY) with respect to Y such that SY = 0. 
With G defined as an l — 2 vector of Lagrange multi- 
pliers, the constrained problem can be rewritten as the 
unconstrained problem, minimizing (Y - EY) T W (Y — 
EY)+2G T SY with respect to Y and G. The solutions 
of this minimization problem (refs. 1 and 8) are 

Y = [E t W£]’* 1 £ t WY - [E T WE}~ l S T G (18) 


and 

G = [S[E T W'£]- 1 S T ] _1 S[E T WE\- l E T WY (19) 

where the superscript —1 denotes a matrix inverse. 

For the user interested in a measure of the error in 
the estimates, the covariance matrix of the constrained 
estimates in terms of Y is readily calculated with the 
matrices used in equations (18) and (19). As shown in 
reference 1, the covariance matrix is 


F = s 2 [£ t W £] -1 

x i - s T [s [e t we] - 1 s T ] _ 1 S [E t WE] ~ 1 

( 20 ) 

In equation (20), s 2 is the estimated variance of the 
measurement error defined by (see refs. 1 and 8) 


s 2 = 


1 


n — 2l 


(Y - £Y) t W(Y - EY) (21) 



Figure 1. Illustration of maximum deviation of rational 
spline from line segments joining knots. 


Tension Adjustment Algorithm 


The Frost and Kinzel tension adjustment algorithm 
is given below; detailed equations are given in ap- 
pendix B. The tension adjustment algorithm for inter- 
polation begins by setting all tension parameters P ^ to 
zero and calculating the cubic-spline interpolating func- 
tion. Each interval is individually tested according to 
its own criterion to determine if the tension must be 
adjusted; if additional tension is required, then Pk is 
incremented by one. For each interval, the criterion to 
be met is that the maximum deviation of the ratio- 
nal spline from the line segment joining the knots not 
exceed a user-specified proportion of the length of that 
line segment. (See fig. 1.) If the maximum deviation ex- 
ceeds the user-specified value, the tension is increased. 
Once all the intervals have been tested and appropriate 
tension parameters incremented, then one of two events 
occurs. If none of the tension parameters changed, then 
the procedure is ended; if at least one tension parame- 
ter changed, then a new rational spline is calculated and 
the maximum deviations are tested. The process of cal- 
culating the interpolating rational spline, testing maxi- 
mum deviations, and adjusting tension parameters can 
be iterated until either no parameter adjustment occurs 
or a specified number of iterations have been made. 
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Rational-Spline Approximation Procedure 

The constrained weighted least-squares problem is 
straightforward to solve and requires no iteration. How- 
ever, because the tension adjustment algorithm is iter- 
ative, combining this algorithm with the least-squares 
solution yields an approximation which must be iter- 
ated and is completely analogous to the rational-spline 
interpolation algorithm. 

The approximation procedure consists of the follow- 
ing steps: 

1. Initially set all tension parameters P It to zero. 

2. With equations (18) and (19), solve the constrained 
weighted least-squares problem for the function ( y *) 
and the second derivative ( y £') values at the knots. 

3. Test the maximum deviation of the rational spline 
in each knot interval; if the maximum deviation ex- 
ceeds the corresponding user-specified bound, incre- 
ment the tension parameter by one. 

4. If the maximum number of iterations has been tried 
or if none of the tension parameters have been 
adjusted, then stop. Otherwise return to step 2 for 
the next iteration. 

The initial pass through step 2 with tension param- 
eter values of zero produces the cubic-spline approxi- 
mation (ref. 8) to the data. Later iterations gradually 
tighten the tension to meet the user’s criteria. An ex- 
ample is presented in the following section to illustrate 
the approximation procedure. 

The rational-spline approximation algorithm does 
have two basic limitations. First, in order for the 
least-squares problem to be well-defined, there must be 
at least three data points on each interval defined by 
consecutive pairs of knots. Second, since the tension 
increases by one during each iteration, the algorithm 
may require an excessive number of iterations to meet 
the user-specified criterion of a small deviation. For 
these cases, in which an essentially linear fit is desired, 
the user is advised to examine the results of several 
different large tension values. 

Example 

The data chosen to illustrate the rational-spline 
approximation consist of 84 values of terrain elevations 
taken during a land survey. The measurements were 
equally spaced at 100-ft intervals and were measured 
with an accuracy of 0.1 ft. 

Nine spline knots were located at the two endpoints 
of the independent- variable range and at seven interme- 
diate points. The intermediate knots were located at 
points which appeared to give a reasonable cubic-spline 
fit. The resulting fit is illustrated in figure 2. The fit 
is quite poor to the right of the fifth knot (located at 
x = 25); the standard deviation of the fit is 1.43 ft. 



Figure 2. Cubic spline fit to terrain elevation data. 



Figure 3. Rational spline fit to terrain elevation data. 


One approach to improving the fit is to add knots 
and to fit a new cubic spline. However, for each knot 
added, two additional unknowns (function value and 
second derivative) must be estimated. The fit can be 
improved without any additional knots by adjusting 
tension between the knots. For the terrain elevation 
data, it appeared that adding tension to the intervals to 
the right of the first, fifth, sixth, and eighth knots could 
improve the poor fit over those intervals. For the first 
interval, a fixed tension of 10 was manually selected. 
For the three remaining intervals, the adjustment algo- 
rithm was applied using the allowed deviations from a 
line given in table I. The results in table I and figure 3 
indicate that tension on only the sixth and eighth in- 
tervals was required to give a much-improved fit. The 
standard error of this rational-spline fit was 0.55 ft and 
required eight iterations. 
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TABLE I. RATIONAL SPLINE FIT TO 
TERRAIN ELEVATION DATA 

[Eight iterations required] 



Knot 

abscissa 

Tension 


Tension 

1 

1.0 

10 


Fixed 

2 

4.0 

0 


Fixed 

3 

7.0 

0 


Fixed 

4 

12.0 

0 


Fixed 

5 

25.0 

0 

15 

Adjusted 

6 

49.0 

7 

10 

Adjusted 

7 

68.0 

0 


Fixed 

8 

74.0 

3 

20 

Adjusted 

9 

84.0 





Concluding Remarks 

An algorithm for obtaining a least-squares approx- 
imation to data with a rational spline has been pre- 
sented. The rational spline combines the advantages 
of a cubic function having continuous first and second 
derivatives with the advantages of a function having in- 
dependently variable tension between consecutive pairs 
of knots. Application of the method of least squares to 
the rational spline leads to a flexible, smooth represen- 
tation of experimental data. 


The example presented illustrates the improved fit 
that can be obtained by approximating with a rational 
spline rather than with a cubic spline. Also demon- 
strated in this example are the choices a user has for 
each consecutive pair of knots: to apply no tension, to 
apply fixed tension, or to determine tension with a ten- 
sion adjustment algorithm. Rational- spline approxima- 
tion also requires the user’s judgment in the selection of 
the number of knots, the knot abscissas, and the allowed 
maximum deviations from line segments. The selection 
of these quantities depends on the actual data and on 
the requirements of a particular application. 

The rational-spline approximation algorithm does 
have two basic limitations. First, in order for the 
least-squares problem to be well-defined, there must be 
at least three data points on each interval defined by 
consecutive pairs of knots. Second, since the tension 
increases by one during each iteration, the algorithm 
may require an excessive number of iterations to meet 
the user-specified criterion of a small deviation. For 
these cases, in which an essentially linear fit is desired, 
the user is advised to examine the results of several 
different large tension values. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
July 11, 1984 
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Appendix A 

Matrix Definitions 

The nonzero entries of the matrices E and S used 
in the constrained least-squares solution are presented 
in this appendix. The matrix definitions are based on 
equations (12) and (15). 

Let x k < Xi < x fc +i for i = 1,2, ...,n. Thus, the 
nonzero entries in row i of matrix E are 

Pi,2k-1 = u 
^i,2fc+l = t 

£ «‘ +2 = Ht “ ') 

where u = (x fc+1 - x)/dx k and t = 1 - u. These 
entries form the submatrix of E containing cofactors 
for subinterval k. 


The nonzero entries of the 1 — 2x1 matrix S are 
the cofactors in equation (15). Hence, these entries are 
given by 


$k,2k-l 


1 

dx k 


q H k 

&k,2k = 

dx k 


Sk,2k+1 — I" 

dx k 


1 

dx k + 1 


Sk,2k+2 


(P k + 2 )H k {Pk + 1 + 2 )H k +i 

dx k dx k - |_i 


S kf 2k+3 — 


1 

dx k + 1 


«Sfc,2fc+4 

where k = 1,2,...,/ — 2. 


dXA; + 1 
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Appendix B 

Equations for Tension Adjustment Algo- 
rithm 

The derivation and description of the automatic ad- 
justment algorithm are given in this appendix. The 
description originally given in reference 6 is based on 
the model given in equation (1). Although the algo- 
rithm presented herein is based on the rational-spline 
definition given in equation (12), the same principle is 
applied: find the maximum deviation of the rational 
spline from the line joining the knots. 

Consider the knots located at (xk,yk) and at 
(^fc+ij yfc+i); the objective is to find the maximum per- 
pendicular distance e/t from the rational spline to the 
line joining the knots. (See sketch.) The distance be- 
tween the rational spline and the line is maximized 
at the point (x*,y*) on the rational spline where the 
derivative of the rational spline is parallel to the line. 
The slope of the line is given by 



, = (l/fe+i ~ Vk) = dyk 

(^fc+l %k) dXfc 


(Bl) 


The line and the derivative of the rational spline are 
parallel if the slopes are equal, that is, if 

F ' k {x) = g (B2) 

Substituting equation (13) into equation (B2) and 


rearranging yields 

\3u 2 {P k t + l)+u 3 P k 


— Vk ~ H k 


+ H k 


(. Pkt + iy 

3 1 2 {P k u + 1) + t 3 P k 


(Pku+iy 

We can substitute u = 1 — t to obtain 

[3(1 -<) 2 (. P k t +!) + (! — tfPk 


Vk + 2/fc+i 

2/fc+i — dyk = o 


-Vk-H k 

+ y k + 1 + H k 

-dy k = 0 


{Pkt + 1 r 

3t 2 [P fc ( l-t) + l] + t 3 P fc 

[Pk{ 1 - t) + l] 2 


- 1 


y'l 


mi 

- 1 j vl 


2/fc+i 

(B3) 


Equation (B3) can be solved for t(0 < t < 1) by any 
one of several iterative methods. One of the simplest 
methods, the secant method (ref. 9), converges to the 
appropriate zero (£*) in just a few iterations. 

After t* is found, the definition of t is used to relate 
x * to t* as follows: 

t . _ {jp - a*) 

dxk 


or 

x* = Xk + t*dxk 

Thus, the ordinate of the point of maximum deviation 
is given by 

y* = F k (x*) 

Now, let Mi = y' k and L\ be the slope and length, 
respectively, of the line joining the knots; let M 2 = 
(y* - yk)/(x* - Xk) and L 2 be the slope and length, 
respectively, of the line joining (xk,yk) and (x*,y*). 
(See sketch.) Since M\ = tan#i and M 2 = tan 0 2 , the 
angle between the two lines 9 = 0 2 — 6\ can be found 
from 

tan# = tan(#2 — #1) 

tan #2 — tan 9 1 
1 + tan #1 tan 9 2 
M 2 — Mi 
1 + Mi M 2 

Finally, the maximum deviation is given by 


e& = L 2 sin 9 


where 

X y/ 2 

L 2 = [(x* - ifc) 2 + (y* - yk ) 2 j 

To use the algorithm, let the user-specified criterion 
be a percent of the length L\. Then, if lOOe^/Li is 
less than or equal to the user-specified percent, Pk does 
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APPENDIX B 


not need to be adjusted; however, if lOOe^/Li is larger algorithm is required. The least-squares algorithm is 

than the user-specified percent, then P % is incremented iterated until user-specified criteria are met on all the 

by one and an additional iteration of the least-squares knot intervals. 
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Symbols 

Ak, Bk, Ck,Dk coefficients of rational spline on 
interval k 

dxk difference between abscissas of 

knots defining interval k 

dyi c difference between ordinates of 

knots defining interval k 

E matrix of cofactors in rational 

spline 

eic maximum deviation of rational 

spline from line for interval k 

E(x) rational spline over entire data set 

F k {x) rational spline on interval k 

G vector of Lagrange multipliers 

I identity matrix 

Li,L 2 lengths of lines 

/ number of knots 

Mi , M 2 slopes of lines 

n number of data points 

Pk tension parameter of rational 

spline on interval k 

S matrix of cofactors in constraint 

equation 

s 2 estimated variance of measurement 

error 


t,u variables used to define rational 

spline 

V covariance matrix of estimates 

W diagonal matrix of weights 


w weight 

x independent variable 

zth abscissa of data 
Xk abscissa of fcth knot 

Y column vector containing ordinates 
of data 

Y column vector containing ordinates 
of data and second derivatives at 
knots 

y dependent variable 

yi 2th ordinate of data 

yk ordinate of kth knot 

0, @1,02 angles used in calculating maxi- 

mum deviation of rational spline 
from line 

Subscripts: 

i zth data point 

k quantity at A;th knot 

Superscripts: 

T matrix transpose 

— 1 matrix inverse 


quantity at point on rational spline 
having maximum deviation from 
line 

A bar over a quantity denotes that quantity at a 
knot. A prime indicates first derivative with respect to 
the independent variable. A double prime indicates sec- 
ond derivative with respect to the independent variable. 
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